stanEx1-1
options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)execute mcmc sampling
goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
mcmc=mdl$sample(
data=data,
seed=1,
chains=4,
iter_sampling=smp,
iter_warmup=wrm,
thin=th,
refresh=1000
)
mcmc
}see mcmc result and parameters
seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
print(mcmc)
prs=mcmc$metadata()$model_params[-1] # reject lp__
for(pr in prs){
if(grepl('^y',pr)) next # not show predictive value "y*" information
if(exc!='' && grepl(paste0('^',exc),pr)) next
drw=mcmc$draws(pr)
if(ch){
par(mfrow=c(1,4))
drw[,1,] |> plot(type='l',main='chain1',ylab=pr)
drw[,2,] |> plot(type='l',main='chain2',ylab=pr)
drw[,3,] |> plot(type='l',main='chain3',ylab=pr)
drw[,4,] |> plot(type='l',main='chain4',ylab=pr)
}
par(mfrow=c(1,2))
drw |> hist(main=pr,xlab='')
drw |> density() |> plot(main=pr)
}
}distributions
normal dist.
distribution of mean, it's from central limit theorem
y~N(m,s), y(-Inf,Inf), E[y]=m, V[y]=s^2ex1.stan
data{
int N;
vector[N] y;
}
parameters{
real m;
real<lower=0> s;
}
model{
y~normal(m,s);
}y=rnorm(10,2,1)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.41 0.71 1.97 1.95 2.95 3.96par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=length(y),y=y)
mdl=cmdstan_model('./ex1.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -95.0671
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 11 -6.8958 4.18942e-05 6.12098e-05 1 1 12
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.3 seconds.mle## variable estimate
## lp__ -6.90
## m 1.95
## s 1.21mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.mcmc## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -7.86 -7.50 1.22 0.85 -10.27 -6.75 1.00 695 778
## m 1.93 1.94 0.53 0.46 1.07 2.75 1.00 643 658
## s 1.51 1.41 0.46 0.37 0.97 2.36 1.00 1219 916mcmc$metadata()$stan_variables## [1] "lp__" "m" "s"mcmc$metadata()$model_params## [1] "lp__" "m" "s"seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -7.86 -7.50 1.22 0.85 -10.27 -6.75 1.00 695 778
## m 1.93 1.94 0.53 0.46 1.07 2.75 1.00 643 658
## s 1.51 1.41 0.46 0.37 0.97 2.36 1.00 1219 916
Bernoulli dist.
The event with probablity p occur y=1, not occur y=0
y~Ber(p), y{0,1}, E[y]=p, V[y]=p(1-p)ex0-0.stan
data{
int N;
array[N] int y;
}
parameters{
real<lower=0,upper=1> p;
}
model{
p~beta(1,1);
y~bernoulli(p);
}data=list(N=9,y=sample(0:1,9,replace=T))
mdl=cmdstan_model('./ex0-0.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -8.08344
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 5 -4.76736 0.00254713 5.10058e-05 1 1 8
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -4.77
## p 0.78mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.mcmc## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -6.94 -6.68 0.68 0.32 -8.34 -6.45 1.00 1018 1579
## p 0.73 0.74 0.13 0.12 0.49 0.91 1.00 870 1055mcmc$metadata()$stan_variables## [1] "lp__" "p"mcmc$metadata()$model_params## [1] "lp__" "p"seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -6.94 -6.68 0.68 0.32 -8.34 -6.45 1.00 1018 1579
## p 0.73 0.74 0.13 0.12 0.49 0.91 1.00 870 1055
Poisson dist.
The event occur l times in unit range, the event occur y times.
y~Po(l), y{0,1...Inf}, E[y]=l, V[y]=lex2-2.stan
data{
int N;
array[N] int y;
}
parameters{
real<lower=0> l;
}
model{
y~poisson(l);
}
generated quantities{
array[N] int y1;
for(i in 1:N)
y1[i]=poisson_rng(l);
}n=10
y=rpois(n,3)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.00 2.25 3.00 3.10 3.75 6.00par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=n,y=y)
mdl=cmdstan_model('./ex2-2.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -1.70601
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 5 4.07347 0.000452794 8.34729e-05 1 1 7
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ 4.07
## l 3.10
## y1[1] 2.00
## y1[2] 2.00
## y1[3] 3.00
## y1[4] 4.00
## y1[5] 5.00
## y1[6] 6.00
## y1[7] 3.00
## y1[8] 2.00
## y1[9] 4.00
## y1[10] 3.00mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ 4.71 4.98 0.68 0.33 3.30 5.22 1.00 1065 1456
## l 3.18 3.14 0.57 0.58 2.31 4.17 1.00 599 1172
## y1[1] 3.21 3.00 1.82 1.48 1.00 7.00 1.00 1716 1732
## y1[2] 3.13 3.00 1.91 1.48 0.00 7.00 1.00 1848 1644
## y1[3] 3.18 3.00 1.85 1.48 1.00 7.00 1.00 1958 1843
## y1[4] 3.11 3.00 1.82 1.48 0.00 6.00 1.00 1687 2091
## y1[5] 3.17 3.00 1.91 1.48 0.00 7.00 1.00 1756 1788
## y1[6] 3.24 3.00 1.90 1.48 1.00 7.00 1.00 1691 1804
## y1[7] 3.16 3.00 1.86 1.48 0.00 6.00 1.00 1900 1819
## y1[8] 3.07 3.00 1.84 1.48 0.00 6.00 1.00 1715 1728
## y1[9] 3.19 3.00 1.89 1.48 1.00 7.00 1.00 1623 1852
## y1[10] 3.07 3.00 1.82 1.48 1.00 6.00 1.00 2003 2047
y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 
binomial dist.
The event with probablity p occur y times on n trials
y~B(n,p), y{0,1...n}, E[y]=np, V[y]=np(1-p)ex2-3.stan
data{
int N;
array[N] int n;
array[N] int y;
}
parameters{
real<lower=0> p;
}
model{
y~binomial(n,p);
}n=10
n0=floor(runif(n,1,10))
y=rbinom(n,n0,0.3)
data=list(N=n,n=n0,y=y)
mdl=cmdstan_model('./ex2-3.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -48.5806
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 9 -21.9007 0.00100885 0.000178818 1 1 10
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -21.90
## p 0.26mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -23.68 -23.42 0.65 0.29 -25.08 -23.20 1.01 882 853
## p 0.28 0.27 0.07 0.07 0.17 0.39 1.01 715 811
log normal dist.
logalithm of variable follows normal distribution
log y~N(m,s), y>0ex2-4.stan
data{
int N;
vector[N]<lower=0> y;
}
parameters{
real m;
real<lower=0> s;
}
model{
y~lognormal(m,s);
}n=10
y=rlnorm(n,2,1)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3.0 6.2 14.8 20.4 22.6 72.1par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=n,y=y)
mdl=cmdstan_model('./ex2-4.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -27.1969
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 12 -13.6408 0.000112385 5.38307e-05 0.9702 0.9702 20
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -13.64
## m 2.58
## s 0.95mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.mcmc## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -14.83 -14.46 1.20 0.84 -16.97 -13.73 1.01 581 694
## m 2.57 2.57 0.41 0.36 1.89 3.24 1.01 698 630
## s 1.18 1.12 0.36 0.28 0.77 1.81 1.00 1211 1139mcmc$metadata()$stan_variables## [1] "lp__" "m" "s"mcmc$metadata()$model_params## [1] "lp__" "m" "s"seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -14.83 -14.46 1.20 0.84 -16.97 -13.73 1.01 581 694
## m 2.57 2.57 0.41 0.36 1.89 3.24 1.01 698 630
## s 1.18 1.12 0.36 0.28 0.77 1.81 1.00 1211 1139
exponential dist.
The event occur l times in unit range, the event take time y to occur
y~ex(l), y>0, E[y]=1/l, V[y]=1/l^2ex2-5.stan
data{
int N;
vector[N] y;
}
parameters{
real<lower=0> l;
}
model{
y~exponential(l);
}n=10
y=rexp(n,1)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.14 0.44 0.84 1.37 1.32 4.82par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=n,y=y)
mdl=cmdstan_model('./ex2-5.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -13.2728
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 4 -13.1256 9.0604e-05 6.08284e-07 1 1 7
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -13.13
## l 0.73mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.mcmc## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -13.90 -13.63 0.72 0.33 -15.33 -13.39 1.00 897 1159
## l 0.81 0.78 0.24 0.23 0.45 1.27 1.00 588 709mcmc$metadata()$stan_variables## [1] "lp__" "l"mcmc$metadata()$model_params## [1] "lp__" "l"seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -13.90 -13.63 0.72 0.33 -15.33 -13.39 1.00 897 1159
## l 0.81 0.78 0.24 0.23 0.45 1.27 1.00 588 709
gamma dist.
The event occur l times in unit range, the event take time y to occur a times
y~Ga(a,l), y>0, E[y]=a/l, V[y]=a/l^2ex2-6.stan
data{
int N;
vector[N] y;
}
parameters{
real<lower=0> a;
real<lower=0> l;
}
model{
y~gamma(a,l);
}n=10
y=rgamma(n,2,1)
summary(y)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.34 0.54 1.01 1.59 2.09 4.72par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')
data=list(N=n,y=y)
mdl=cmdstan_model('./ex2-6.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -25.7033
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 8 -14.1624 0.00023761 3.35903e-05 1 1 11
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -14.16
## a 1.53
## l 0.96mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.mcmc## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -14.47 -14.13 1.07 0.73 -16.54 -13.50 1.00 558 902
## a 1.98 1.89 0.73 0.71 0.95 3.29 1.00 461 501
## l 1.31 1.24 0.54 0.52 0.56 2.35 1.01 464 427mcmc$metadata()$stan_variables## [1] "lp__" "a" "l"mcmc$metadata()$model_params## [1] "lp__" "a" "l"seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -14.47 -14.13 1.07 0.73 -16.54 -13.50 1.00 558 902
## a 1.98 1.89 0.73 0.71 0.95 3.29 1.00 461 501
## l 1.31 1.24 0.54 0.52 0.56 2.35 1.01 464 427
beta dist.
use as prior of binomial dist parameter p
p~Be(a,b), p[0,1], E[p]=a/(a+b), V[p]=ab/(a+b)^2/(a+b+1)ex2-7.stan
data{
int N;
vector[N] p;
}
parameters{
real<lower=0> a;
real<lower=0> b;
}
model{
p~beta(a,b);
}n=20
p=rbeta(n,1,2)
summary(p)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.002 0.173 0.300 0.368 0.560 0.871par(mfrow=c(1,2))
hist(p,main='p')
density(p) |> plot(main='p')
data=list(N=n,p=p)
mdl=cmdstan_model('./ex2-7.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -6.81721
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 8 2.50745 6.66508e-05 0.000171474 0.8139 0.8139 12
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ 2.51
## a 0.89
## b 1.56mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.mcmc## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ 1.91 2.27 1.11 0.74 -0.30 2.92 1.00 821 872
## a 0.98 0.96 0.27 0.26 0.58 1.44 1.01 741 1030
## b 1.75 1.71 0.52 0.50 0.99 2.70 1.01 692 869mcmc$metadata()$stan_variables## [1] "lp__" "a" "b"mcmc$metadata()$model_params## [1] "lp__" "a" "b"seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ 1.91 2.27 1.11 0.74 -0.30 2.92 1.00 821 872
## a 0.98 0.96 0.27 0.26 0.58 1.44 1.01 741 1030
## b 1.75 1.71 0.52 0.50 0.99 2.70 1.01 692 869
dirichlet dist.
use as prior of categorical dist parameter p
p~dir(p0), p[0,1]ex2-8.stan
data {
int N;
int K;
matrix[N, K] p;
}
parameters {
simplex[K] p0;
}
model {
for(i in 1:N){
p[i]~dirichlet(p0);
}
}library(gtools)
n=20
p=rdirichlet(n,c(0.5,0.3,0.2))
summary(p)## V1 V2 V3
## Min. :0.002 Min. :0.000 Min. :0.000
## 1st Qu.:0.054 1st Qu.:0.020 1st Qu.:0.006
## Median :0.358 Median :0.128 Median :0.055
## Mean :0.404 Mean :0.292 Mean :0.304
## 3rd Qu.:0.724 3rd Qu.:0.501 3rd Qu.:0.539
## Max. :0.999 Max. :0.998 Max. :0.951boxplot(p)
data=list(N=n,K=ncol(p),p=p)
mdl=cmdstan_model('./ex2-8.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = 53.6938
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 6 60.8296 0.000136795 3.11352e-07 1 1 9
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ 60.83
## p0[1] 0.45
## p0[2] 0.32
## p0[3] 0.24mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.mcmc## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ 56.44 56.76 1.07 0.69 54.39 57.38 1.00 997 1240
## p0[1] 0.44 0.44 0.06 0.06 0.35 0.53 1.00 2254 1496
## p0[2] 0.32 0.32 0.06 0.05 0.23 0.41 1.01 1909 1036
## p0[3] 0.24 0.24 0.05 0.05 0.17 0.32 1.00 2184 1495mcmc$metadata()$stan_variables## [1] "lp__" "p0"mcmc$metadata()$model_params## [1] "lp__" "p0[1]" "p0[2]" "p0[3]"seeMCMC(mcmc)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ 56.44 56.76 1.07 0.69 54.39 57.38 1.00 997 1240
## p0[1] 0.44 0.44 0.06 0.06 0.35 0.53 1.00 2254 1496
## p0[2] 0.32 0.32 0.06 0.05 0.23 0.41 1.01 1909 1036
## p0[3] 0.24 0.24 0.05 0.05 0.17 0.32 1.00 2184 1495
